#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import glob, sys
import numpy as np
import xarray as xr
import pandas as pd

## FUEL CONSTANTS ##
# fuelheat   = 8000.    # Combustion heat of drey fuel (Btu/lb)
fueldens   = 32.      # Ovendry particle density (lb/ft^3)
st         = 0.0555   # Fuel total mineral content (1)
se         = 0.010    # Fuel effective mineral content (1)

def ros_calculation(fmc_ij,GSI,U,f,fuelmce_live):
    ## FUEL ARRAY ##
    # Fuel models and fuel sizes [20,6]
    w0_ij = np.array([[0.10,  0.00,  0.00,  0.00,  0.00, 1.00,  0.00], #GR2
                      [0.25,  0.00,  0.00,  0.00,  0.00, 1.90,  0.00], #GR4
                      [0.10,  0.00,  0.00,  0.00,  0.00, 3.40,  0.00], #GR6
                      [1.00,  0.00,  0.00,  0.00,  0.00, 5.40,  0.00], #GR7
                      [0.50,  1.00,  0.00,  0.00,  0.00, 7.30,  0.00], #GR8
                      [1.00,  1.00,  0.00,  0.00,  0.00, 9.00,  0.00], #GR9
                      [1.35,  2.40,  0.75,  0.00,  0.00, 0.00,  3.85], #SH2
                      [0.45,  3.00,  0.00,  0.00,  0.00, 0.00,  6.20], #SH3
                      [3.60,  2.10,  0.00,  0.00,  0.00, 0.00,  2.90], #SH5
                      [3.50,  5.30,  2.20,  0.00,  0.00, 0.00,  3.40], #SH7
                      [2.05,  3.40,  0.85,  0.00,  0.00, 0.00,  4.35], #SH8
                      [4.50,  2.45,  0.00,  0.00,  0.00, 1.55,  7.00], #SH9
                      [0.20,  0.90,  1.50,  0.00,  0.00, 0.20,  0.90], #TU1
                      [0.95,  1.80,  1.25,  0.00,  0.00, 0.00,  0.20], #TU2
                      [1.10,  0.15,  0.25,  0.00,  0.00, 0.65,  1.10], #TU3
                      [4.00,  4.00,  3.00,  0.00,  0.00, 0.00,  3.00], #TU5
                      [0.50,  2.20,  2.80,  0.00,  0.00, 0.00,  0.00]])*0.0459137  #TL3 Fuel load (lb/ft)

    savr_ij = np.array([[2000.,  109.,  30.,  8.,  1800., 1800.,  1500.], #GR2
                        [2000.,  109.,  30.,  8.,  1800., 1800.,  1500.], #GR4
                        [2200.,  109.,  30.,  8.,  2000., 2000.,  1500.], #GR6
                        [2000.,  109.,  30.,  8.,  1800., 1800.,  1500.], #GR7
                        [1500.,  109.,  35.,  8.,  1300., 1300.,  1500.], #GR8
                        [1800.,  109.,  30.,  8.,  1600., 1600.,  1500.], #GR9
                        [2000.,  109.,  30.,  8.,  2000., 2000.,  1600.], #SH2
                        [1600.,  109.,  30.,  8.,  2000., 2000.,  1400.], #SH3
                        [750.0,  109.,  30.,  8.,  2000., 2000.,  1600.], #SH5
                        [750.0,  109.,  30.,  8.,  2000., 2000.,  1600.], #SH7
                        [750.0,  109.,  30.,  8.,  2000., 2000.,  1600.], #SH8
                        [750.0,  109.,  30.,  8.,  1800., 1800.,  1500.], #SH9
                        [2000.,  109.,  30.,  8.,  1800., 1800.,  1600.], #TU1
                        [2000.,  109.,  30.,  8.,  2000., 2000.,  1600.], #TU2
                        [1800.,  109.,  30.,  8.,  1600., 1600.,  1400.], #TU3
                        [1500.,  109.,  30.,  8.,  2000., 2000.,  750.0], #TU5
                        [2000.,  109.,  30.,  8.,  2000., 2000.,  1500.]])#TL3  Surface-area-to-volume ratio (1/ft)

    # Fuel models [20]
    fuelheat = np.array([8.0, 8.0, 8.0, 8.0, 8.0,        # GR2, GR4, GR6, GR7, GR8
                         8.0, 8.0, 8.0, 8.0, 8.0,        # GR9, SH2, SH3, SH5, SH7
                         8.0, 8.0, 8.0, 8.0, 8.0,        # SH8, SH9, TU1, TU2, TU3
                         8.0, 8.0])*1000.                # TU5, TL3 Combustion heat of dry fuel (Btu/lb)

    fueldepth = np.array([1.00, 2.00, 1.50, 3.00, 4.00,  # GR2, GR4, GR6, GR7, GR8
                          5.00, 1.00, 2.40, 6.00, 6.00,  # GR9, SH2, SH3, SH5, SH7
                          3.00, 4.40, 0.60, 1.00, 1.30,  # SH8, SH9, TU1, TU2, TU3
                          1.00, 0.20])                   # TU5, TL3 Fuel depth (ft)

    fuelmce = np.array([15., 15., 40., 15., 30.,       # GR2, GR4, GR6, GR7, GR8
                        40., 15., 40., 15., 15.,       # GR9, SH2, SH3, SH5, SH7
                        40., 40., 20., 30., 30.,       # SH8, SH9, TU1, TU2, TU3
                        25., 30.])                     # TU5, TL3 Moisture content of extinction (%)

    waf = np.array([0.36, 0.42, 0.39, 0.46, 0.49,      # GR2, GR4, GR6, GR7, GR8
                    0.52, 0.36, 0.44, 0.55, 0.55,      # GR9, SH2, SH3, SH5, SH7
                    0.46, 0.50, 0.33, 0.36, 0.38,      # SH8, SH9, TU1, TU2, TU3
                    0.36, 0.28])                       # TU5, TL3 Wind adjustment factor  

    # Fuel sizes [6]
    #savr_ij = np.array([2000., 109., 30., 8., 2000., 1500.])                 # Surface-area-to-volume ratio (1/ft)
    dead = [0,1,2,3,4]
    live = [5,6]

    ## CURING ##
    T = 1./(1.-0.5)-GSI/(1.-0.5)
    if T>1.:
        T=1.
    w0_ij[f,4]=T*w0_ij[f,5]
    w0_ij[f,5]=w0_ij[f,5]-T*w0_ij[f,5]

    ## WEGHTING FACTORS ##
    A_ij = savr_ij[f,:]*w0_ij[f,:]/fueldens
    A_dead = np.sum(A_ij[dead])
    A_live = np.sum(A_ij[live])
    A_T = np.sum(A_ij)
    if (A_dead>0.) & (A_live>0.):
        f_ij = np.concatenate([A_ij[dead]/A_dead, A_ij[live]/A_live])
    elif (A_dead==0.):
        f_ij = np.concatenate([A_ij[dead], A_ij[live]/A_live])
    elif (A_live==0.):
        f_ij = np.concatenate([A_ij[dead]/A_dead, A_ij[live]])
    f_dead = A_dead/A_T
    f_live = A_live/A_T

    dead_bool = np.array([True, True, True, True, True, False, False])
    live_bool = np.array([False, False, False, False, False, True, True])
    g_ij = f_ij.copy()
    g_ij[(savr_ij[f,:]>=1200.) & (dead_bool)] = np.sum(f_ij[dead][savr_ij[f,dead]>=1200.])
    g_ij[(savr_ij[f,:]>=1200.) & (live_bool)] = np.sum(f_ij[live][savr_ij[f,live]>=1200.])
    g_ij[(savr_ij[f,:]<1200.) & (savr_ij[f,:]>=192.) & (dead_bool)] = np.sum(f_ij[dead][(savr_ij[f,dead]<1200.) & (savr_ij[f,dead]>=192.)])
    g_ij[(savr_ij[f,:]<1200.) & (savr_ij[f,:]>=192.) & (live_bool)] = np.sum(f_ij[live][(savr_ij[f,live]<1200.) & (savr_ij[f,live]>=192.)])
    g_ij[(savr_ij[f,:]<192.) & (savr_ij[f,:]>=96.) & (dead_bool)] = np.sum(f_ij[dead][(savr_ij[f,dead]<192.) & (savr_ij[f,dead]>=96.)])
    g_ij[(savr_ij[f,:]<192.) & (savr_ij[f,:]>=96.) & (live_bool)] = np.sum(f_ij[live][(savr_ij[f,live]<192.) & (savr_ij[f,live]>=96.)])
    g_ij[(savr_ij[f,:]<96.) & (savr_ij[f,:]>=48.) & (dead_bool)] = np.sum(f_ij[dead][(savr_ij[f,dead]<96.) & (savr_ij[f,dead]>=48.)])
    g_ij[(savr_ij[f,:]<96.) & (savr_ij[f,:]>=48.) & (live_bool)] = np.sum(f_ij[live][(savr_ij[f,live]<96.) & (savr_ij[f,live]>=48.)])
    g_ij[(savr_ij[f,:]<48.) & (savr_ij[f,:]>=16.) & (dead_bool)] = np.sum(f_ij[dead][(savr_ij[f,dead]<48.) & (savr_ij[f,dead]>=16.)])
    g_ij[(savr_ij[f,:]<48.) & (savr_ij[f,:]>=16.) & (live_bool)] = np.sum(f_ij[live][(savr_ij[f,live]<48.) & (savr_ij[f,live]>=16.)])
    g_ij[savr_ij[f,:]<16.]=0.
    g_ij[savr_ij[f,:]<16.]=0.

    ## MOISTURE OF EXTINCTION ##
    if (A_live!=0.):
        W = np.sum(w0_ij[f,dead]*np.exp(-138./savr_ij[f,dead]))/np.sum(w0_ij[f,live]*np.exp(-500./savr_ij[f,live]))
        Mfdead = np.sum(fmc_ij[dead]*w0_ij[f,dead]*np.exp(-138./savr_ij[f,dead]))/np.sum(w0_ij[f,dead]*np.exp(-138./savr_ij[f,dead]))
        fuelmce_live = (2.9*W*(1.-Mfdead/fuelmce[f])-0.226)*100. 
        if fuelmce_live<fuelmce[f]:
            fuelmce_live=fuelmce[f].copy()

    ## CHARACTERISTIC VALUES ##
    wn_ij = w0_ij[f,:]*(1.-st)                        # Net fuel load (lb/ft)
    wn_dead = np.sum(g_ij[dead]*wn_ij[dead])          # Net dead fuel load (lb/ft)
    wn_live = np.sum(g_ij[live]*wn_ij[live])          # Net live fuel load (lb/ft)
    Mf_dead = np.sum(f_ij[dead]*fmc_ij[dead])         # Dead fuel moisture content
    Mf_live = np.sum(f_ij[live]*fmc_ij[live])         # Live fuel moisture content
    savr_dead = np.sum(f_ij[dead]*savr_ij[f,dead])      # Dead fuel surface-area-to-volume ratio (1/ft)
    savr_live = np.sum(f_ij[live]*savr_ij[f,live])      # Live fuel surface-area-to-volume ratio (1/ft)
    savr = f_dead*savr_dead+f_live*savr_live          # Characteristic surface-area-to-volume ratio (1/ft)
    bulkdens = np.sum(w0_ij[f,:])/fueldepth[f]        # Mean bulk density (lb/ft^3)
    beta = bulkdens/fueldens                          # Mean packing ratio (1)
    beta_op = 3.348*savr**(-0.8189)                   # Optimum packing ratio (1)
    ns = 0.174*se**(-0.19)                                                                                    # Mineral damping coefficient 
    nm_dead = 1.-2.59*(Mf_dead/fuelmce[f])+5.11*(Mf_dead/fuelmce[f])**2.-3.52*(Mf_dead/fuelmce[f])**3.        # Dead fuel moisture damping coefficient
    nm_live = 1.-2.59*(Mf_live/fuelmce_live)+5.11*(Mf_live/fuelmce_live)**2.-3.52*(Mf_live/fuelmce_live)**3.  # Live fuel moisture damping coefficient
    if nm_dead<0.:
        nm_dead=0.
    if nm_live<0.:
        nm_live=0.

    ## HEAT SOURCE ##
    gamma_max = (savr**1.5)/(495.+0.0594*(savr**1.5))                    # Maximum reaction velocity (1/min)
    A = 133.*savr**(-0.7913)                                             # Coef for optimum reaction velocity
    gamma = gamma_max*np.exp(A*(1.-(beta/beta_op)))*(beta/beta_op)**A    # Optimum reaction velocity (1/min)
    Ir = gamma*fuelheat[f]*ns*(wn_dead*nm_dead+wn_live*nm_live)          # Reaction intensity (btu/ft^2 min)
    flux = np.exp((0.792+0.681*savr**0.5)*(beta+0.1))/(192.+0.2595*savr) # Propagating flux ratio
    heatsource = Ir*flux

    ## HEAT SINK ## 
    Qig_ij = 250.+11.16*fmc_ij     # Heat of preignition (btu/lb)
    heatsink = bulkdens*(f_dead*np.sum(f_ij[dead]*np.exp(-138./savr_ij[f,dead])*Qig_ij[dead])+f_live*np.sum(f_ij[live]*np.exp(-138./savr_ij[f,live])*Qig_ij[live]))

    ## WIND FACTOR ##
    C = 7.47*np.exp(-0.133*savr**0.55)
    B = 0.02526*savr**0.54
    E = 0.715*np.exp(-3.59e-4*savr)
    U_adj = U*waf[f]*196.850
    if U_adj>(0.9*Ir):               # Wind limit
        U_adj=0.9*Ir
    phiW = C*(U_adj)**B*(beta/beta_op)**(-E)  

    ## SPREAD RATE ##
    ros0 = heatsource/heatsink  
    ros = ros0*(1.+phiW)
    return ros, Mf_dead, Mf_live, fuelmce_live

def main(index,date,lat,lon,model,fuel,contribution,anom_prec):

    index_str=str(index).zfill(4)
    year=date.split('-')[0]
    month=date.split('-')[1]
    day=date.split('-')[2]
    fuelmce_live = np.array([15., 15., 40., 15., 30.,       # GR2, GR4, GR6, GR7, GR8
                             40., 15., 40., 15., 15.,       # GR9, SH2, SH3, SH5, SH7
                             40., 40., 20., 30., 30.,       # SH8, SH9, TU1, TU2, TU3
                             25., 30.])                     # TU5, TL3 Moisture content of extinction (%)

    ## PRESENT ##
    path_fmc='../OUTS/FMC_Nelson_csv/'
    df_fmc_1h = pd.read_csv(path_fmc+'N'+index_str+'_fmc1h.csv',parse_dates=['time'])
    df_fmc_1h = df_fmc_1h.loc[(df_fmc_1h.time>=str(date+' 12:00:00')) & (df_fmc_1h.time<=str(date+' 20:00:00'))]
    df_fmc_10h = pd.read_csv(path_fmc+'N'+index_str+'_fmc10h.csv',parse_dates=['time'])
    df_fmc_10h = df_fmc_10h.loc[(df_fmc_10h.time>=str(date+' 12:00:00')) & (df_fmc_10h.time<=str(date+' 20:00:00'))]
    df_fmc_100h = pd.read_csv(path_fmc+'N'+index_str+'_fmc100h.csv',parse_dates=['time'])
    df_fmc_100h = df_fmc_100h.loc[(df_fmc_100h.time>=str(date+' 12:00:00')) & (df_fmc_100h.time<=str(date+' 20:00:00'))]
    path_lfmc='../OUTS/LFMC/ERA5/'
    ds_fmc_herb = xr.open_dataarray(path_lfmc+'Herb/ERA5_fmch_'+year+'-'+month+'.nc').sel(time=date)      
    ds_fmc_woody = xr.open_dataarray(path_lfmc+'Woody/ERA5_fmcw_'+year+'-'+month+'.nc').sel(time=date)
    path_uv='../DATA/ERA5/'
    ds_u10 = xr.open_dataarray(path_uv+'U10m/ERA5_10m_u_component_of_wind_'+year+'-'+month+'.nc').sel(time=date)      
    ds_v10 = xr.open_dataarray(path_uv+'V10m/ERA5_10m_v_component_of_wind_'+year+'-'+month+'.nc').sel(time=date)
    path_gsi='../OUTS/GSI/ERA5/means/'
    ds_gsi = xr.open_dataset(path_gsi+'ERA5_gsi_'+year+'-'+month+'.nc').GSI.sel(time=date)      

    fmc_1h = df_fmc_1h['FMC'].mean()*100.
    fmc_10h = df_fmc_10h['FMC'].mean()*100.
    fmc_100h = df_fmc_100h['FMC'].mean()*100.
    fmc_1000h = df_fmc_100h['FMC'].mean()*100.
    fmc_herb = ds_fmc_herb.interp(latitude=lat, longitude=lon).squeeze()
    fmc_woody = ds_fmc_woody.interp(latitude=lat, longitude=lon).squeeze()
    fmc = np.array([fmc_1h,
                    fmc_10h,
                    fmc_100h,
                    fmc_1000h,
                    fmc_1h,
                    fmc_herb.data,
                    fmc_woody.data])
    u10 = ds_u10.interp(latitude=lat, longitude=lon).squeeze()
    v10 = ds_v10.interp(latitude=lat, longitude=lon).squeeze()
    U=np.sqrt(u10**2.+v10**2.).mean()
    GSI = ds_gsi.interp(latitude=lat, longitude=lon).squeeze()
    del fmc_1h, fmc_10h, fmc_100h, fmc_1000h, fmc_herb, fmc_woody

    ros_f1, fmd_f1, fml_f1, fmel_f1 = ros_calculation(fmc,GSI.data,U,fuel,fuelmce_live[fuel],'control')

    ## PICONTROL ##
    path_lfmc='../OUTS/LFMC/ERA5/'
    ds_fmc_herb = xr.open_dataarray(path_lfmc+'Herb/ERA5_fmch_'+year+'-'+month+'.nc').sel(time=date)      
    ds_fmc_woody = xr.open_dataarray(path_lfmc+'Woody/ERA5_fmcw_'+year+'-'+month+'.nc').sel(time=date)
    path_gsi='../OUTS/GSI/ERA5/means/'
    ds_gsi = xr.open_dataset(path_gsi+'ERA5_gsi_'+year+'-'+month+'.nc').GSI.sel(time=date)      

    if contribution=='all':
        fmc_1h = df_fmc_1h['FMC_piControl_'+model].mean()*100.
        fmc_10h = df_fmc_10h['FMC_piControl_'+model].mean()*100.
        fmc_100h = df_fmc_100h['FMC_piControl_'+model].mean()*100.
        fmc_1000h = df_fmc_100h['FMC_piControl_'+model].mean()*100.
        del ds_fmc_herb, ds_fmc_woody, ds_gsi
        path_lfmc='../OUTS/LFMC/'+model+'/piControl/'
        ds_fmc_herb = xr.open_dataarray(path_lfmc+'Herb/ERA5_fmch_'+year+'-'+month+'.nc').sel(time=date)      
        ds_fmc_woody = xr.open_dataarray(path_lfmc+'Woody/ERA5_fmcw_'+year+'-'+month+'.nc').sel(time=date)
        path_gsi='../OUTS/GSI/'+model+'/piControl/means/'
        ds_gsi = xr.open_dataset(path_gsi+'ERA5_gsi_'+year+'-'+month+'.nc').GSI.sel(time=date)      
    elif contribution=='dead':
        fmc_1h = df_fmc_1h['FMC_piControl_'+model].mean()*100.
        fmc_10h = df_fmc_10h['FMC_piControl_'+model].mean()*100.
        fmc_100h = df_fmc_100h['FMC_piControl_'+model].mean()*100.
        fmc_1000h = df_fmc_100h['FMC_piControl_'+model].mean()*100.
    elif contribution=='live':
        fmc_1h = df_fmc_1h['FMC'].mean()*100.
        fmc_10h = df_fmc_10h['FMC'].mean()*100.
        fmc_100h = df_fmc_100h['FMC'].mean()*100.
        fmc_1000h = df_fmc_100h['FMC'].mean()*100.
        del ds_fmc_herb, ds_fmc_woody, ds_gsi
        path_lfmc='../OUTS/LFMC/'+model+'/piControl/'
        ds_fmc_herb = xr.open_dataarray(path_lfmc+'Herb/ERA5_fmch_'+year+'-'+month+'.nc').sel(time=date)      
        ds_fmc_woody = xr.open_dataarray(path_lfmc+'Woody/ERA5_fmcw_'+year+'-'+month+'.nc').sel(time=date)
        path_gsi='../OUTS/GSI/'+model+'/piControl/means/'
        ds_gsi = xr.open_dataset(path_gsi+'ERA5_gsi_'+year+'-'+month+'.nc').GSI.sel(time=date)      


    fmc_herb = ds_fmc_herb.interp(latitude=lat, longitude=lon).squeeze()
    fmc_woody = ds_fmc_woody.interp(latitude=lat, longitude=lon).squeeze()
    fmc_piControl = np.array([fmc_1h,
                    fmc_10h,
                    fmc_100h,
                    fmc_1000h,
                    fmc_1h,
                    fmc_herb.data,
                    fmc_woody.data])
    GSI_piControl = ds_gsi.interp(latitude=lat, longitude=lon).squeeze()
    del fmc_1h, fmc_10h, fmc_100h, fmc_1000h, fmc_herb, fmc_woody

    ros_piControl_f1, fmd_piControl_f1, fml_piControl_f1, fmel_piControl_f1 = ros_calculation(fmc_piControl,GSI_piControl.data,U,fuel,fmel_f1,contribution)

    ## SSP245 ##
    path_lfmc='../OUTS/LFMC/ERA5/'
    ds_fmc_herb = xr.open_dataarray(path_lfmc+'Herb/ERA5_fmch_'+year+'-'+month+'.nc').sel(time=date)      
    ds_fmc_woody = xr.open_dataarray(path_lfmc+'Woody/ERA5_fmcw_'+year+'-'+month+'.nc').sel(time=date)
    path_gsi='../OUTS/GSI/ERA5/means/'
    ds_gsi = xr.open_dataset(path_gsi+'ERA5_gsi_'+year+'-'+month+'.nc').GSI.sel(time=date)      

    if contribution=='all':
        fmc_1h = df_fmc_1h['FMC_ssp245_'+model].mean()*100.
        fmc_10h = df_fmc_10h['FMC_ssp245_'+model].mean()*100.
        fmc_100h = df_fmc_100h['FMC_ssp245_'+model].mean()*100.
        fmc_1000h = df_fmc_100h['FMC_ssp245_'+model].mean()*100.
        del ds_fmc_herb, ds_fmc_woody, ds_gsi
        path_lfmc='../OUTS/LFMC/'+model+'/ssp245/'
        ds_fmc_herb = xr.open_dataarray(path_lfmc+'Herb/ERA5_fmch_'+year+'-'+month+'.nc').sel(time=date)      
        ds_fmc_woody = xr.open_dataarray(path_lfmc+'Woody/ERA5_fmcw_'+year+'-'+month+'.nc').sel(time=date)
        path_gsi='../OUTS/GSI/'+model+'/ssp245/means/'
        ds_gsi = xr.open_dataset(path_gsi+'ERA5_gsi_'+year+'-'+month+'.nc').GSI.sel(time=date)      
    elif contribution=='dead':
        fmc_1h = df_fmc_1h['FMC_ssp245_'+model].mean()*100.
        fmc_10h = df_fmc_10h['FMC_ssp245_'+model].mean()*100.
        fmc_100h = df_fmc_100h['FMC_ssp245_'+model].mean()*100.
        fmc_1000h = df_fmc_100h['FMC_ssp245_'+model].mean()*100.
    elif contribution=='live':
        fmc_1h = df_fmc_1h['FMC'].mean()*100.
        fmc_10h = df_fmc_10h['FMC'].mean()*100.
        fmc_100h = df_fmc_100h['FMC'].mean()*100.
        fmc_1000h = df_fmc_100h['FMC'].mean()*100.
        del ds_fmc_herb, ds_fmc_woody, ds_gsi   
        path_lfmc='../OUTS/LFMC/'+model+'/ssp245/'
        ds_fmc_herb = xr.open_dataarray(path_lfmc+'Herb/ERA5_fmch_'+year+'-'+month+'.nc').sel(time=date)      
        ds_fmc_woody = xr.open_dataarray(path_lfmc+'Woody/ERA5_fmcw_'+year+'-'+month+'.nc').sel(time=date)
        path_gsi='../OUTS/GSI/'+model+'/ssp245/means/'
        ds_gsi = xr.open_dataset(path_gsi+'ERA5_gsi_'+year+'-'+month+'.nc').GSI.sel(time=date)      

    fmc_herb = ds_fmc_herb.interp(latitude=lat, longitude=lon).squeeze()
    fmc_woody = ds_fmc_woody.interp(latitude=lat, longitude=lon).squeeze()
    fmc_ssp245 = np.array([fmc_1h,
                    fmc_10h,
                    fmc_100h,
                    fmc_1000h,
                    fmc_1h,
                    fmc_herb.data,
                    fmc_woody.data])
    GSI_ssp245 = ds_gsi.interp(latitude=lat, longitude=lon).squeeze()

    del ds_fmc_herb, ds_fmc_woody, ds_gsi, fmc_1h, fmc_10h, fmc_100h, fmc_1000h, fmc_herb, fmc_woody

    ros_ssp245_f1, fmd_ssp245_f1, fml_ssp245_f1, fmel_ssp245_f1 = ros_calculation(fmc_ssp245,GSI_ssp245.data,U,fuel,fmel_f1,contribution)

    return ros_f1, fmd_f1, fml_f1, fmel_f1, ros_piControl_f1, fmd_piControl_f1, fml_piControl_f1, fmel_piControl_f1, ros_ssp245_f1, fmd_ssp245_f1, fml_ssp245_f1, fmel_ssp245_f1    

if __name__ == '__main__':
    main(int(sys.argv[1]), str(sys.argv[2]) ,float(sys.argv[3]), float(sys.argv[4]), str(sys.argv[5]), int(sys.argv[6]), str(sys.argv[7]))
